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The Kuramoto model for an ensemble of coupled oscillators provides a paradigmatic example 
of non-equilibrium transitions between an incoherent and a synchronized state. Here we analyze 
populations of almost identical oscillators in arbitrary interaction networks. Our aim is to extract 
topological features of the connectivity pattern from purely dynamical measures, based on the fact 
that in a heterogeneous network the global dynamics is not only affected by the distribution of the 
natural frequencies, but also by the location of the different values. In order to perform a quantitative 
study we focused on a very simple frequency distribution considering that all the frequencies are 
equal but one, that of the pacemaker node. We then analyze the dynamical behavior of the system 
at the transition point and slightly above it, as well as very far from the critical point, when it is 
in a highly incoherent state. The gathered topological information ranges from local features, such 
as the single node connectivity, to the hierarchical structure of functional clusters, and even to the 
entire adjacency matrix. 

PACS number(s): 89.75.-k, 89. 75. He, 05.45.Xt 



I. INTRODUCTION 

Nowadays, it is widely acknowledged that complex pat- 
terns of interaction are ubiquitous in nature as in so- 
ciety Nonetheless, further research is required to 
completely understand how the topology affects the sys- 
tem dynamics [H, Q . In particular how global dynamical 
properties are related with the units dynamics and the 
interactions between them. A unique answer cannot be 
provided since complex networks respond differently de- 
pending on the dynamical processes that take place on 
them 

One of the most interesting of these macroscopi- 
cally defined dynamical processes is synchronization, an 
emerging phenomenon in which populations of interact- 
ing units display a common periodic behavior [l, Q . In- 
deed, understanding the role of connectivity in synchro- 
nization has been the subject of intense research in recent 
years Q • On the one hand, much work has focused on the 
generic properties of dynamical systems, mainly looking 
for necessary and sufficient conditions that would grant 
that a population of units under a set of simple dynam- 
ical rules is able to synchronize Q. On the other hand, 
much progress has been made by studying precise models 
of phase oscillators, being one of the most paradigmatic 
the model proposed by Kuramoto [loj , where the in- 
teraction between the units is proportional to the sine of 
the phase difference. 

In the present work, we will continue along this line 
and analyze a population of Kuramoto oscillators with a 
precise distribution of frequencies. The original work by 
Kuramoto and many subsequent studies considered that 
the oscillators, each coupled equally to all the others, had 
natural frequencies taken from a given distribution. The 
non-zero width of those distributions made the units fol- 
low different trajectories, whereas the interaction term 
made their phases approach. In fact and depending on 



the width of the frequency distribution, there is a critical 
value of the interaction strength above which the units 
tend to entrain their phases and hence leave the incoher- 
ent regime. If the natural frequencies of the oscillators 
are identical, a unique outcome is possible as the only 
attractor of the dynamics is a completely synchronized 
state in which all the oscillators end up in a common 
phase. And this occurs for any initial conditions and for 
any (connected) topology [36j . 

In systems with regular patterns of connectivity (in- 
cluding all-to-all) the only complexity comes from the 
frequency distribution, whereas in more realistic (non- 
homogeneous) patterns, not only the frequency values 
matter but the precise location as well [1, [ll| . 

Here we will focus on a particular frequency distribu- 
tion, one which is just one step away from the homoge- 
neous case. Such distribution has identical frequencies for 
all oscillators except one. This singular oscillator, with 
a higher frequency than the rest, has received the name 
of pacemaker and its effect in populations of Kuramoto 
oscillators has been analyzed [l3. Il3| . In [l^l, Kori and 
Mikhailov consider a special case where the pacemaker 
affects its neighbors but it is not affected by them; under 
these conditions they find numerically that the range of 
frequencies of the pacemaker for which the system can 
attain global synchronization depends on the " depth" of 
the network; defining the depth as the maximum dis- 
tance from the pacemaker to peripheral nodes. Radicchi 
and Meyer- Ortmanns [Tsj consider regular structures for 
which the conditions to synchronize can be analytically 
computed. 

In this paper we use several properties of the hetero- 
geneity induced by the existence of the pacemaker to find 
useful relations between topology and dynamics. On one 
hand, by knowing the topology one should be able to 
infer the dynamical properties of the network. On the 
other hand, by measuring the dynamics some structural 
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properties can be inferred, and this will be our purpose. 

First, we use a similar procedure than the one used 
in and [l^l, showing that there is a critical value 
for the frequency of the pacemaker above which the (fre- 
quency) synchronized state cannot exist. This is related 
with the existence of a synchronized solution (also ex- 
ploited in that applies to any subset of oscillators. 
Wc find, however, that from a practical point of view the 
most restrictive condition is usually for the equation of 
the pacemaker that involves its connectivity, and hence 
there is a clear relationship between the critical frequency 
and the pacemaker connectivity which can be used as an 
experimental measure of the degree. 

In order to get more details on the network structure 
we analyze the system above the critical value where cor- 
relations between dynamical evolution of the nodes ap- 
pear. Such correlations enable to reveal the hierarchical 
organization and to recover the network connectivity. 

The structure of the paper is as follows. First, in Sec. 2, 
we characterize the coherent state and the transition to 
the incoherent one by means of a proper definition of the 
order parameter. Then, in Sec. 3, we qualitatively an- 
alyze the behavior of the system when it is not in the 
frequency-locked state. Sec. 4 is devoted to study the re- 
lation between local connectivity and the ability of the 
system to reach a synchronized (frequency-locked) state. 
In Sec. 5 we focus on the system slightly above the tran- 
sition towards the incoherent state. We show that it is 
possible to perform some hierarchical analysis concerning 
the connectivity network. Finally, in Sec. 6 we study the 
system far above the critical point, in a regime character- 
ized by short range correlations where it becomes easy to 
identify the nodes directly connected to the pacemaker. 
Thus the reconstruction of the whole connectivity pat- 
tern is accurate and fast. 



II. SYNCHRONIZATION AND PHASE 
TRANSITION 

In the original Kuramoto model [1, , the phases of 
the oscillators evolve according to the following equation 

N 

(fii = oJi + (J^sm{(pj - Lpi), (1) 

where N is the total number of units of the system, uji 
is the natural frequency of unit i, taken from a distri- 
bution, and tr stands for the coupling strength. This 
case corresponds to a fully connected topology, i.e. each 
unit interacts with all the other ones. The ability of the 
system to reach a coherent state, for a given coupling 
strength, depends only on the width of the distribution 
of natural frequencies. 

Here we want to consider arbitrary connectivity pat- 
terns. In this situation, the behavior of the system can 
no longer be understood in terms of the ratio between 
the distribution width and the coupling strength only. It 



is also relevant where the natural frequencies values are 
located, since on a generic interaction network nodes are 
not equivalent anymore. 

From now on we are using the 2 levels hierarchical net- 
work of 9 nodes represented in Fig.[T]as a benchmark and, 
when not otherwise stated, all the figures refer to that 
connection pattern. This network has been presented 
in [l5| as a very simple example of the class of determin- 
istic scale-free hierarchical networks proposed by Ravasz 
and Barabasi in [l^ . We choose this small regular con- 
nectivity pattern as a simple paradigmatic example show- 
ing general properties of the studied systems, since it 
makes easy to recognize the role of each node. 




FIG. 1: (Color online) Hierarchic network that will be used 
as benchmark. In this particular setting the pacemaker is 
located on a peripheral node (marked as P) of degree kp =3. 
The other nodes are grouped into sets using different colors. 
The elements of each set are topologically equivalent if we 
look at the network from the point of view of the pacemaker. 
Consequently their dynamical evolution is identical. 

Let us rewrite the equation for the evolution of the 
phases including a connectivity matrix Oy- that is sym- 
metric and takes values 1(0) if node i and j are connected 
(disconnected): 

N 

ipi = uji + ^ Qij sm{ipj ~ ipi), (2) 

where we have rescaled time by setting a — I. Now 
wc consider all the oscillators to have the same natural 
frequency (0 without loss of generality), except one of 
them, called the pacemaker, whose frequency is w 7^ 0. 
It is precisely this extremely simple choice of frequen- 
cies that enables to study the roles played by individual 
oscillators. 

If a stationary state exists, then all the effective fre- 
quencies will take constant values and the following con- 
ditions have to be satisfied: 

JV 

aij sm{(pj - ifi) = fli Vi ^ p (3) 

i=i 

N 

CO + ^ Up J sm{ipj — ipp) = flp (4) 

where {fii} are the effective frequencies of the oscillators. 
Notice that summing up eqs. ©-(HI) the coupling terms 
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cancel because of the symmetry of the interaction and it 
results in 

JV 

J2 = ^- (5) 

i=l 

Looking at eqs. ©-(HI) it is easy to recognize that there 
is an interplay between two effects. On the one hand the 
width of the frequencies distribution (in our present case 
this role is played by w itself) tends to keep the evolution 
of the oscillators apart since each one follows its natu- 
ral frequency. On the other hand, the interaction term 
makes them to approach their phases as well as their 
effective frequencies. Then we conclude that if the pace- 
maker natural frequency is small enough, the interaction 
term dominates and, after a transient time, all effective 
frequencies 17,; will be identical 

n, = Lo/N \fi, (6) 

including the pacemaker. In this case we can say that 
the system is in a frequency-locked state, since all oscil- 
lators have the same frequency although the phases are 
not equal, because there is a coupling term (that of the 
pacemaker) that cannot vanish. 

When increasing the pacemaker frequency w, some 
oscillators cannot keep the phase difference and the 
frequency-locked state is broken. The left-hand side 
of eq. is indeed bounded because of the sine terms, 
whereas the right term increases as the pacemaker fre- 
quency is increased. A similar conclusion can be deduced 
from eq. Q. Consequently, there will be a transition 
from a synchronized to an incoherent state. Thus we 
can define the critical value Wp as the maximum value 
of the natural frequency of the pacemaker for which the 
system can attain global synchronization. 

Such a transition for a population of phase oscillators is 
typically characterized by an order parameter R, defined 
through the equation: R e** = J^j ^^'^^ i where 4' is a 
global phase (not constant) [l7| . 

Anyway, in the present work, following [ll|, [l^, we 
adopt another order parameter that is a normalized mea- 
sure of the effective frequency dispersion (standard devi- 
ation) : 

where (w) is the average effective frequency of the oscil- 
lators population, a constant quantity always equal to 
Lo/N. According to its definition, takes values in the 
interval [0 , 1] (see Fig. [5]). It should be noticed that, since 
above the critical frequency the system is not able to 
reach a steady state anymore, calculation of the order 
parameter ([7]) requires to perform averages over an ap- 
propriate time window. Anyway, the value of (w) does 
not change because what we found in ([5|) is a general 
result, even for instantaneous values of the effective fre- 
quencies. 



To find the precise value of the critical frequency we 
apply the Newton-Raphson method (NR) and check, as 
a function of the frequency a;, whether the synchronized 
solution of eqs. ©-(HI) exists. To simulate the dynam- 
ics of the system in the incoherent state (w > w^) we 
take as initial phases {vj(0)} the stationary values of 
the differences provided by the NR solution for w = oj^. 
Eqs. ^ are numerically integrated with Eulcr's Method 
(first-order), unless otherwise stated, at fixed time step 
5t ^ 10-2. 
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FIG. 2: (Color online) Order parameter as a function 
of the natural frequency of the pacemaker. Different curves 
correspond to different settings: (•) refers to the pacemaker 
located on nodel in Fig.[l] (A;p = 8), (a) to the pacemaker 
on node 2 (fcp = 3). The average value (ri^)* for uj > ujp was 
calculated on a time window At = 100. 



III. INCOHERENT STATE 

Above the critical frequency the system is no longer 
in a stationary state and hence the effective frequencies 
arc no longer constant. 

Numerical simulations show that, after a transient 
time, the system enters into a "periodic" state (see Fig[3|). 
The features of this periodic state are not affected by the 
initial conditions and they only depend on the pacemaker 
frequency and location. It is precisely this fact that en- 
ables to infer topological properties from dynamical mea- 
surements. 

Fig. 2] summarizes what we have learned up to now, 
shedding light on some interesting details. The time av- 
erage of the effective frequency of the pacemaker {(pp)t 
and that of one of its neighbor {^Pj)t are plotted as func- 
tions of the pacemaker natural frequency. These quan- 
tities are calculated from numerical simulations taking 
into account appropriate time windows. 

Starting from small values of w, the picture shows how 
all the effective frequencies increase together linearly, fol- 
lowing the reference line Vti = lo/N defined by eq. ([5]). 
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FIG. 3: (Color online) Effective frequencies above the crit- 
ical point as functions of time. That of the pacemaker 
(red, top curve), in this particular setting located on node 3 
in Fig.[T] is on average much larger than the others (lower 
curves). Panels a, b, c and d refer respectively to a pacemaker 
natural frequency value that is 1.01, 1.05, 1.2 and 2 times its 
critical value. Time starts after a transient lag Ta — 10. 



Then, when uj reaches the critical value w^, they do sep- 
arate. Initially, the average effective frequency of the 
pacemaker goes through a more than linear increasing, 
while the others start decreasing, keeping their (average) 
values very close to each other. For even larger values, 
when u! ^ ujp , Fig. |4] shows how the average effective fre- 
quency {'fp)t tends to w, asymptotically increasing along 
a new reference line with slope equal to 1. At the same 
time, {(pi}t for i ^ p goes to zero, as required by the 
conservation law (El). 



IV. CRITICAL FREQUENCY AND LOCAL 
TOPOLOGY 



In this section we explore the relation between the 
topology of the network and the value of the critical nat- 
ural frequency of the pacemaker depending on the node 
where it is located. 

Let us begin by writing the equation for the pacemaker 
in the synchronized state. As a consequence of eq. ([S]), we 
have 



N 



J = l 



, sin((/3j — ipp) = uj/N. 



(8) 



This equation links the natural frequency of the pace- 
maker to the constant values of the phase differences be- 
tween it and its neighbors, when all the units are oscillat- 
ing with the same effective frequency. Since the number 
of non-null terms ajp in the previous expression is given 
by the number of nodes connected with the pacemaker 
and s'm(fj — fi) £ [—1, 1], the degree (or connectivity) 




FIG. 4; (Color online) Average effective frequencies as func- 
tion of the natural frequency of the pacemaker uj. The figure 
shows the behavior of two oscillators: node 1 (black circles) 
and node 2 (red triangles). On the right side the pacemaker 
is node 1, on the left one it is node 2. Initially the frequencies 
are synchronized and they increase linearly with slope 1/N 
(dashed line) as expected from eq. Then, when lj reaches 
the critical value, that is different for different location of the 
pacemaker, they get apart. Far above the critical values, the 
average frequency of the pacemakers approach asymptotically 
a new reference line with slope 1 (solid line). The time aver- 
ages were performed on a time window At — 100. 



of the pacemaker is a bound for the absolute value of the 
sum in eq. 

Thus there is an upper bound for the critical frequency: 



N 



(9) 



where kp is the degree of the pacemaker. Indeed, any 
value larger than the right term in the inequality ^ is 
surely unable to satisfy eq. ([5]) and hence the system is 
unable to be frequency synchronized. 

Notice that we have obtained this bound by taking into 
account a single equation, that of the pacemaker. We can 
write for any oscillator the analogous of eq. ^ as follows 



sin(iy9j — ipi) ~ (^/N, Vi ^ p. 



(10) 



It is easy to verify that no stricter condition can arise 
from any of these equations jSTj . However, stronger 
bounds could exist due to the combination of eq. ([5]) and 
some of eqs. pH)) . 

Let us consider a set of (n + l) connected nodes, among 
which the pacemaker is included [s^ • Labeling them by 
an increasing index z = l,2,...,n + l= p and summing 
up their equations we obtain: 



(n + l) 



UJ 

N 



n+l N 



2^ 2^ fly sin(v3j 



(11) 
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If two nodes in the considered group are neighbors their 
respective interaction terms cancel each other. So the 
number of remaining terms of the sums in cq. (fTTj) is given 
by: 



n+l 



n+l 

E 



(12) 



where ki is the degree of the th node and Kout is equal 
to the number of links connecting the nodes of the con- 
sidered set to external ones. 

Consequently, cq. (|lip can be rewritten as: 



(n + l)- =0;- 



E sin(0/), 



(13) 



1=1 



where = Lpj — Lpi, being i and j connected nodes re- 
spectively inside and outside the group. 

We are now able to write the expression of the upper 
bound for the critical frequency in a generalized form: 



< < Kout ^ 



N 



1 



= N 



(14) 



where Nout stands for the number of nodes not belonging 
to the considered set. Eq. ([H)) reduces to the previous 
upper bound if one chooses n = 0. 

In this way we can write a very large number of condi- 
tions, that is the number of the connected sets of nodes 
that include the pacemaker and which size ranges from 
1 to iV — 1. Among these, the strongest one is that for 
which the ratio Kout I Nout takes its minimum value. This 
is a combinatorial problem, in principle very simple, but 
hard from a computational point of view, since the num- 
ber of conditions grows at least exponentially with the 
network size. 

Minimizing the ratio Kout /Nout we find the strictest 
condition on that can be expressed in the form of a 
single equation. No other equation obtained as a linear 
combination of equations ©-(HI) may provide a stronger 
bound. This condition is analogous to the necessary con- 
dition for global synchronization concerning the surface 
(here Kout) of any subset of nodes derived in IJ] for 
randomly distributed natural frequencies and generic os- 
cillators. However, these conditions are not sufficient. 
In our case, it is not sure that the Kout remaining sine 
terms of eq. (|13p are allowed to take their minimal values 
simultaneously. This kind of problems directly involves 
the sine functions arguments that may be not indepen- 
dent since they are differences between pairs of phases 
and we are dealing with a system of N coupled equa- 
tions. It may happen that two or more phases are tied 
among each other by a certain set of equations of the 
kind fi(ipi,{ipi.}) = (where the nodes {ij} are neigh- 
bors of the node i). Consequently we cannot minimize 
the sum of sine terms on a hypercube [0; 2tt] °"* but we 
have to restrict ourselves on a hyper-surface of dimen- 
sion Kout — "m, where m is the number of constrains. A 



system may experience this kind of difficulty (that we 
can regard as a kind of angles frustration) only if cycles 
are present and there is some anisotropy, and only when 
1 < kp < N — 1. Therefore, for a good number of regu- 
lar connectivity patterns, as those analyzed in [l3| . there 
is not such a problem and it is possible to analytically 



calculate the all set of values {wi^'}, p = 1 

As a simple analytically solvable network let us con- 
sider a Caylcy tree with coordination number z, made 
up of S shells. For each node it is indeed possible to sin- 
gle out a connected "group" such that Kout = 1, taking 
in all the nodes on the branch starting from the consid- 
ered pacemaker. In this way we are minimizing the ratio 
Kout /Nout so that we can consider the strictest equation 
among cqs. ([Ti]) . Moreover, since there is not any cycle, 
neither there are problems of angles frustration. There- 
fore, the obtained expressions give the correct values, not 
just bounds. In this way we obtain for the critical fre- 
quency: 



= N 



1 



N-Etoi^-^y 



where s is the shell of the pacemaker. 

Even though in real complex networks it is not so easy 
to calculate the {uj^}-, we have empirically verified that 
only in few cases the critical frequency is much smaller 
than its first upper bound This can be clearly ob- 
served in Fig. [51 where we plotted the ratios between the 
real critical values and the corresponding upper bound, 
for every choice of the pacemaker in several networks. 

The accuracy of this estimation enables us to use it in 
the opposite direction, i.e. to get an estimation of the 
pacemaker degree from an experimental measure of the 
critical frequency. We can invert eq. ([9]) obtaining: 



kr) ^ 



,iV- 1 



(15) 



but, since the right term is not an integer, the smallest 



allowed value for kp is 



k: 



,N-1 
' N 



(16) 



where [x] stands for the integer part of x. We can con- 
clude that eq. ([T6)) gives the correct value of kp whenever 



(kp - 1) 



N 



k 



N 



iV-1' ^7V-1 



This fact implies that the estimator for the degree of 
the pacemaker is very reliable. Indeed, it only fails when 
the critical frequency is really smaller than its bound ([9|) . 



V. SLIGHTLY ABOVE THE CRITICAL POINT 

In this and in the next section we translate the rich 
dynamical information that the system provides in the 
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FIG. 5: (Color online) Critical frequency of a pacemaker as a 
function of its degree, for a set of networks. We have divided 
the critical frequency by the degree and by N/{N — 1) such 
that the bound given by eq. Q is 1. We have shifted the 
data for the different networks and the horizontal lines are 
the reference (equal to 1) for each case. From bottom to top 
the networks are: 1) Zachary club social network [l^ used 
in community detection applications; 2) hierarchical network 
of 125 nodes and 3 levels [l^; 3) network of 4 communities 
of 32 nodes each used as benchmark in community detection 
algorithms [l^ where all the nodes have the same degree; 
4) network of jazz bands [23|; 5-6-7) three networks of three 
levels of community structure used to relate topological and 
temporal scales in synchronization [2l[; 8) Caenorhabditis 
elegans neural network [2^ . 



incoherent state into useful topological information. Here 
we focus on the behavior of the system slightly above the 
critical point, while in Sec. 6 we will analyze the system 
when the natural frequency of the pacemaker is many 
time larger than its critical value. 

We are interested in estimating how much similar two 
nodes are from a global topological perspective. For 
this purpose we need to define an appropriate correla- 
tion function, able to relate the dynamical responses of 
pairs of oscillators. 

Looking for the expression of a good correlation func- 
tion, we get no help from the average values {ipi)t = 
ipi{t)dt. Indeed, in this regime, all the oscillators, 
except the pacemaker, have the same average effective 
frequency. On the contrary, it can be useful to look at 
the difference between instantaneous values. We mea- 
sure the frequency of every oscillator at each time, inside 
a suitable interval. In order to define a correlation, that 
is a quantity that has to be non-negative and symmetric 
with respect to the two nodes indexes i and j, it is rea- 
sonable to start from a power of the absolute value of the 
difference \(pl^\t) — (p'j'\t)\, where (p) stand for the pace- 
maker that induces the considered dynamical evolution. 



Therefore, we propose 



Dividing by w makes that the argument of the root is less 
than 1 because, even if the frequencies may take nega- 
tive values (see Fig. [3]), the condition \(pi{t)\ <^ uj always 
holds. 

The period of the effective frequencies oscillation de- 
pends on which node is the pacemaker. Then, in order 
to compute averages on time that are really independent 
from the considered interval, we have to choose a time 
window many times larger than the oscillation period. 
Furthermore, since \ipi — Lpp\ ^ \(pi — ipj\ for any i,j ^ p, 
we decide to exclude these contributions, taking into ac- 
count only terms of the kind cf ^ where i ^ p and j ^ p. 

Finally, in order to remove the dependence from the 
index p we have to average on all the possible pacemakers. 
Summarizing in a compact expression, our correlation 
function can be written as follow: 



1 



N 

E 



1 




N-2 ^,..ti-hjt 



A. Hierarchical organization 



■dt 



(17) 



Once we have obtained the correlation matrix we 
can proceed to some hierarchical analysis. In the 
present work we use the standard Unweighted Pair Group 
Method Average (UPGMA) [H algorithm to compute 
such diagrams. What we find out is a hierarchy of dy- 
namical communities, whose meaning is immediately un- 
derstandable in the case of small networks, such as our 
benchmark in Fig. [T] (see Fig. [51). 

Obviously, this simple network docs not need any 
analysis to obtain its hierarchical organization, but this 
methodology can be very useful when applied to func- 
tional hierarchical network. 

As a paradigmatic example, let us consider the cor- 
ticocortical network of the cat at the macroscopic level. 
We look at each cortical area as a basic unit, modeling 
it as a Kuramoto oscillator, finding out similar results as 
in[2i,p. 

In Fig. [7] we show that, going down along our den- 
dogram starting from the root, it is possible to recog- 
nize two communities clearly separated. Then, the right 
branch splits up into two parts and the left one undergoes 
into two subsequent bifurcations, so that it is possible to 
identify three groups of nodes on it. At this level we 
have five communities. Four of them correspond to well 
known physiological sub-systems: the fronto-limbic (FL), 
the somatosensory- motor (SM), the auditory (A) and the 
visual (V). The fifth one (HUBS) is composed - except 
for a single area [s^- by super-hubs, sometimes consid- 
ered as a meta-community (rich-club) [l^, HH . The most 
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FIG. 6: (Color online) The network of Fig.[T] and its corre- 
sponding dendrogram. Correlations are calculated averaging 
on a time window At — 60, after a transient lag Ta = 10, for 

LJ = 1.1 ■ UJp. 




n 

FIG. 7: (Color online) Dendrogram of the cortical brain net- 
work of the cat. Different colors correspond to different sub- 
systems: the fronto- limbic (FL), the somatosensory-motor 
(SM), the auditory (A) and the visual (V). The rich-club is 
labeled with HUBS, while the branch indicated with the label 
hp (pink) is the area that belongs to the hippocampus and it 
is out of its place. Correlations are calculated averaging on 
a time window At = 100, after a transient lag Ts = 10, for 
a; = 1.1 • Wp. 



relevant aspect of our hierarchical analysis is that there 
is no way to recognize this meta-community if the den- 
drogram is constructed by means of static methods. Nei- 
ther it can be obtained throughout correlation matrices 
generated from the adjacency matrix using, for instance, 
Pearson's Coefficient |26l| . Nor these nodes emerge as a 
community when the modularity function is maximized. 
Indeed, maximizing the modularity we obtain as an op- 



timal partition the same 4 groups corresponding to the 4 
physiological sub-systems. 

In general, complex networks can be organized, and 
thus analyzed, at different hierarchical levels. For social 
networks it is very important that a group is tight, so 
that the multiple connections within the group give rise 
to the concept of community. On the contrary, in biolog- 
ical networks the most crucial concept is function rather 
than connectivity per se. Therefore, methods that rely 
on the connections within groups and maximize modular- 
ity will not be enough to identify biological units, based 
primarily on function (27l. [28l|. In this case, our method, 
which analyzes the dynamical correlation between units, 
provides a better approach to infer functional relation- 
ships. 

One of the known problems of the methods commonly 
used for detecting community structures in complex net- 
works is the existence of the so called resolution limit, 
found by Fortunato and Barthelemy [2^. This issue is 
related to the impossibility for the methods based on 
modularity optimization to go beyond certain resolution 
which is related to the community size and to the num- 
ber of links between communities. The paradigmatic ex- 
ample of the problem is a network formed by "cliques" 
(small groups of totally connected nodes) which are very 
sparsely connected among them. We have checked such 
structures and found that dynamically the correlations 
are very strong within the cliques and not among nodes 
belonging to different modules, showing that our method 
detecting the hierarchical organization is not affected by 
the resolution limit problem. 



B. Recovering network topology 

Let us now take a step backward and recover some- 
thing we had previously discarded. In the sum of equa- 
tion (|17p we had excluded terms in which one of the in- 
dexes was equal to p since they were heterogeneous. So 
Cij = jfZ2 '^p^i j '^j ■ Aiiy^^'Yi ^-Iso the set of elements 
cFpp p = 1, . . . , iV contains information. We may ask our- 
selves which are the oscillators most strongly correlated 
with the pacemaker and if they share some topological 
property. The simplest hypothesis is that the set of kp 
largest c^^- identifies the neighbors of the pacemaker. This 
is reasonable since, even if the pacemaker is very weakly 
correlated with the rest of the oscillators, coefficients c^^ 
are not uniform and the topological distance is the most 
immediate quantity we may suppose this variability is re- 
lated to. In Sec 4 we showed how to find out an estimator 
of the degree of each node from the critical frequencies. 
Thus if we are able to select the possible neighbors we 
would be in principle able to reconstruct the entire net- 
work. 

The first problem we face in the attempt to validate 
this hypothesis is that our list of likely neighbors gives 
us an asymmetric and weighted adjacency matrix, which 
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elements are 

a'pj. = for i = fc* + 1, . . . ,iV, 

where k* is the estimator for the degree of the pacemaker 
given by and c^^- . > c^j^ whenever i < k* and I > k*. 

Moreover ^ since generally speaking c™„ 7^ 
c"j„. Therefore we have to remove the weights and to 
symmetrize this matrix. Here we propose an algorithm 
to perform this task that is at the same time simple and 
efficient. It consists in four steps. 

1) Symmetrize the matrix in the usual way: = 

2) Compute a list of temporary degree k'^ > fc* as the 
number of non-null elements a^^; 

3) Order all the non-zero values a^„j in a list, from the 
smaller to the larger; 

4) Check which ones among the corresponding likely 
links have to be removed, starting from the weakest one. 

We proceed as follows: given a pair of nodes m and 
n whose link is the weakest one, if and only if k'^ > k^ 
and fc'j > A:* we remove that link, setting = a^„j = 

0. In this case both fc^ and k'„ are reduce by one unit. 
Otherwise we go to the next link, going on along the 
entire list, till the strongest link. 

This method roots in the hypothesis, empirically very 
well verified, that the matrix a^„ contains all the links 
of the real network, plus a number of false positive ones, 

1. e. that there is no false negative link. Thus we need 
just to remove, never to add edges. 

Moreover it works properly only if our estimators {fc* } 
of the actual degrees {kn} are correct, otherwise we may 
make additional errors. Fortunately it is a very infre- 
quent problem. The sole hypothesis we make is that 
the probability for a link of being a "false" one is a 
monotonously decreasing function of the correlation be- 
tween the nodes it joins. 

Finally, the method does not ensure that in the final 
estimated network fc^ = fc* Vn because it is possible that 
even if k'^ > fc* , the n-th oscillator has no possible neigh- 
bor which temporary degree is larger than its estimated 
one. Sometime this fact may cause new errors, some oth- 
ers it acts as a compensation of the underestimation of 
the real degrees. 

In order to quantify how good a reconstruction is, we 
introduce the following error definition: 

err% = ^P±^ ■ 100, 

where Fp and Fn are respectively the number of false 
positive (spurious) and false negative (missing) links in 
the reconstructed network, and L is the number of edges 
in the original connectivity pattern. Globally speaking, 
we can state that our method allows for a reconstruction 
of an arbitrary connectivity pattern with a good preci- 
sion. Taking into account the networks in TableUl on 
average we have err% = 6.5. 



Among these networks there are artificial as well as real 
connectivity patterns. They were selected to be represen- 
tative of several classes of networks, including hierarchi- 
cal as well as not hierarchical, with and without commu- 
nity structure, regular and not regular. For this reason, 
the average error calculated on this set of benchmarks 
can be considered as a good estimator of the accuracy of 
the proposed reconstruction method when applied on a 
given unknown connectivity pattern. 
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TABLE I: Results of the reconstruction on several networks. 
On the columns we list: the size of the system (N), the total 
number of links in the original network (L), the total error 
in the estimation of the degrees (Kerr = ~ 
total number of links in the reconstructed network before the 
removal of exceeding links (L'), the number of false positive 
(Fp) and false negative (Fn) links in this network, the same 
for the final reduced network (LV, Fp/Fn) and the final total 
error (err%). From the first row, the networks are: our usual 
benchmark [31, ring of 6 cliques of 3 nodes [2^, hierarchical 
network of 25 nodes and 2 levels [l^, Zachary club social 
network [l^, ring of 16 cliques of 3 nodes [2^, cortical brain 
network of the cat [s^], hierarchical network of 125 nodes and 
3 levels [l^, network of 4 communities of 32 nodes each [l^ . 
2 networks of 3 levels of community structure ^J]. 



VI. FAR FROM THE CRITICAL POINT 

Far above the critical point the system behaves quite 
differently. As clearly shown in Fig. [8] all the units are 
characterized by effective frequencies that, after a tran- 
sient time, oscillate around precise values that are equal 
to their own natural frequency, as can be seen in the left 
panels of Fig.|8l From this point of view, by increasing 
the natural frequency of the pacemaker the coupling is 
less and less important. But, in any case, there are still 
reminiscences of the interactions since the amplitudes of 
the oscillations decay very fast with the distance from the 
pacemaker. Indeed, the frequencies of the neighbors of 
the pacemaker oscillate with an amplitude that is roughly 
Aneigh — 2, while all the other oscillators are almost at 
rest if compared with them. These conditions allow us to 
recognize the neighbors of a given pacemaker even if we 
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do not know how many they are. Therefore, we may de- 




Time 



FIG. 8: (Color online) Effective frequencies as function of 
time far above the critical point (lj = 20 ■ lj^) . Plots refer to 
the same network used for the previous pictures, in the case 
of 2 different choice of the pacemaker: node 1 (fci = 8) above; 
node 2 (^2 = 3) below. On the left hand side we plotted the 
frequency of all the nodes in the network. On the right one 
the scale has been changed and the pacemakers are left out. 
Notice how above, where all the nodes are neighbors of the 
pacemaker, we may observe a unique curve. On the contrary, 
below there are 2 different kinds of oscillations. The largest 
ones are those of the neighbors of the pacemaker, the others 
belong to the oscillator not directly connected to it. Time 
starts after a transient lag Ts = I. The integration time step 
used is St — 10"''. 



fine a simplified correlation function that better suits this 
situation and that only connects each pacemaker with its 
neighbors: 



maxt(y)jft)) - iniiit{ipi{t)) _ 
maxt(v3p(t)) - mint((^p(t)) Ap 



(18) 



Previous expression is the ratio between two positive 
terms (amplitudes) and it is equal to 1 for i = p. 

On any connectivity pattern, the values c^^ are dis- 
tributed along a stair whose highest step is easy to iden- 
tify even if we consider short time windows. The tran- 
sient time, indeed, is always very short in this regime. 
We do not need any more to completely reconstruct the 
entire connection topology. 

All we have to do is to compute the values c^j for 
each pacemaker. After finding out the maximum val- 
ues maxi^p Vp, we choose an appropriate threshold, 
say 0.5. A node j will be a neighbor of the pacemaker p 
if / {maxi^p Cp^) > 0.5. Now we are able to construct 
a connectivity matrix. 

Let us notice that in this case there is no need for 
symmetrization since the adjacency matrix constructed 
in this way is already symmetric because this method is 
based on a reliable general property that holds for any 
connectivity pattern. The use of a threshold is therefore 
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Nodes Labels 

FIG. 9: (Color online) Normalized correlations of the cortical 
brain network of the cat for pacemaker on node 1 (fcp = 10). 
The (•) are the correlation values calculated through eq. H17|) 
for uj = 1.2- uj^ (At = 100, = 20). The (a) correspond 
to the correlations given by expression (|18|) when lj = 20 ■ Up, 
calculated on a time window At = 1 and waiting a transient 
time Ts =0.1. All the values have been divided by the maxi- 
mum on each set (excluding the auto-correlation). Notice that 
while in the first case there is an almost continuous spectrum 
of values, in the second one it is easy to identify a group of 
points (in red) above the line at 0.5 clearly separated from 
the rest. Those are the 10 neighbors of node 1. 



in principle unnecessary, since all the neighbors have the 
same amplitude of the frequency oscillation, when the 
pacemaker natural frequency is above a certain value. 
But, since this value is not know a priori and it may 
be very large if the distribution of the degrees among 
the neighbors of the pacemaker is very wide, it is useful 
from an empirical point of view. It is important to stress 
that, even if we are still in a regime where some degree of 
heterogeneity among the neighbors is conserved, there is 
no chance to make any error in the recovered topology. 
Indeed; the amplitudes of the frequency oscillations of 
oscillators not directly connected to the pacemaker are 
at least one order of magnitude smaller than those of its 
neighbors (see Figs. [8] and [9]). By means of this method 
all the topologies considered in Table|T] arc properly re- 
constructed, without errors. 

In addition, not all nodes need to be considered as 
pacemakers. While the method discussed in Sec. V-B 
requires to perform dynamical measures for every possi- 
ble location of the pacemaker, for the current description 
this is not necessary. Indeed, we can look for the neigh- 
bors of a number N' < N of pacemakers in order to 
get all the connections in the considered network. From 
an experimental point of view, adopting the conceptual 
framework proposed in [3]| , we may consider the choice 
of a certain pacemaker as the application of a drift on 
a given unit in a system of identical coupled oscillators. 
This means that it is possible to solve the problem with 
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less than N experiments. 

The criterion for choosing the ordered sequence of 
nodes on which we locate the pacemaker can vary. We 
may operate a random extractions, or we may start from 
a randomly chosen node and then move to one of its 
neighbors along a random walk. Another option, much 
more convenient especially in the case of scale-free net- 
works, can be adopted if the critical frequencies associ- 
ated to each oscillator are known. We can ordered the 
nodes according to decreasing critical frequency, start- 
ing from the highest one. In this way wc proceed from 
larger to smaller (estimated) degrees, taking an impor- 
tant advantage if the degrees distribution is not uniform 
and there are hubs. The hubs, indeed, provide informa- 
tion about a large number of links by means of very few 
experiments fFiglTO|. 



VII. CONCLUSIONS 

Systems of non-identical Kuramoto oscillators have 
been recently shown to display a degree of synchroniza- 
tion that depends strongly on the topology of the under- 
lying complex network. Here, these dynamical proper- 
ties, in particular by setting different types of correlations 
between the dynamical evolution of the oscillators, have 
been used to gather information on the connectivity pat- 
terns. Remarkably, this is the case for most experimen- 
tal situations, where the a priori unknown connectivity 
of a particular network is inferred from purely dynamical 
measurements. 

When the oscillators are identical (all of them having 
the same natural frequency) any topological configura- 
tion has a unique attractor, which is the complete syn- 
chronized state; synchronized meaning that the oscilla- 
tors end up in such a state that all effective frequencies 
and phases are identical. This state does not offer any in- 
formation about the topology. We perturb this setting by 
allowing one of the oscillators to have a different natural 
frequency than the rest. This unit is called the pace- 
maker of the network. Such perturbation causes that 
the final state is no longer phase-synchronized. But if 
the natural frequency of the pacemaker is not very dif- 
ferent from the value of the rest of the population, the 
system still will keep a certain degree of synchronization, 
since the whole system can evolve with the same effective 
frequency. However, if the frequency difference becomes 
larger, the system will be unable to find any kind of syn- 
chronization. The threshold between the former case and 
this latter is a well defined value, which is strictly depen- 
dent on the location of the pacemaker in the network. 
In this context, we can use the correlations between the 
effective frequencies of the oscillators in such incoherent 
state to reproduce the network connectivity. 

Moreover, we show that the dynamical correlations in 
different situations, whether close of far from the critical 
point, provide complementary information on the net- 




# Nodes # Nodes 

FIG. 10: (Color online) Average number of reconstructed 
links as a function of the number of nodes we considered 
as pacemakers (number of trials). From the top to the bot- 
tom, the considered networks are: a pair of Barabasi- Albert 
networks, respectively with parameter A; = 3 (left) and A: = 10 
(right); a pair of Erdos-Reyni graphs with average degree 
equal to 15 (left) and 60 (right); a pair of random regu- 
lar graphs with degree 5 (left) and 100 (right). The size is 
= 1000 for all of them. Different lines corresponds to dif- 
ferent selection algorithms. Blue dashed lines stand for the 
ordered sequence on the basis of the critical frequencies val- 
ues; the red solid ones for the random walk; the green dots 
ones for random extractions. Both the random walk and the 
random extractions are averaged over 1000 samples. The hor- 
izontal black line marks the %90 of links: notice how in any 
case we never need more than %70 of the nodes in order to 
reconstruct %90 of the links, decreasing to a %30 — %40 in 
the case of the ordered sequence for scale-free networks. Cor- 
relations are computed under the same conditions as those of 
Fig.H 



work: 

1. Working around the critical point we are able to 
estimate the degree of each pacemaker merely by 
its critical frequency. 

2. Slightly above the transition point the hierarchical 
structure of the whole network (related to func- 
tional modules) can be obtained from the correla- 
tions between effective frequencies. A further re- 
finement enables to recover the whole connection 
network with a good degree of accuracy. 

3. Far above the critical point it is possible to rec- 
ognize which are the oscillators that are directly 
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connected to an individual pacemaker from a very 
short measurement of the time evolution of the ef- 
fective frequencies. In this way we can recover the 
connectivity pattern and this method turns out to 
be much more precise and more efficient than the 
previous one. 

In summary, this paper deals with different approaches 
relating dynamical properties of individual nodes to the 
topology of the network. The topological properties in- 
ferred from dynamics can be local (the existence of a link 
between two nodes) as well as global (hierarchical organi- 
zation of the nodes in the functional network) . In partic- 
ular, for a scale-free network and if the node degrees are 
known (or have been estimated from the critical frequen- 
cies), considering 30% of the possible pacemakers, always 
selecting the most connected nodes, will be enough to re- 
construct approximately 90% of the links. 

Other papers have considered the reconstruction of the 
network from dynamical information. Similar to our pro- 
posal with specific targets, Tegner et. al. [l^l analyzed 
the dynamical response of a gene-regulatory network by 
changing expression levels of par ticular genes. On the 
contrary, Di Bernardo et. al. [Sj] considered the global 
effect of different types of perturbations to infer the net- 



work topology. This approach has been followed recently 
also by Gorur Shandilya and Timme in [33 |. where it 
is assumed that there is some information about the dy- 
namical evolution of the isolated units and about the cou- 
pling. Our method, based on the change of the frequency 
of a single unit and how it enhances correlations among 
the nodes, can be more effective in oscillatory systems. In 
any case, for practical purposes the method chosen will 
depend on the specific details of the experimental setup 
and even a combination of different ones can be the most 
appropriate. 
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